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ABSTRACT 

We consider height-integrated equations of an advection-dominated accretion flow (ADAF), 
assuming that there is no mass outflow. We include convection through a mixing length for- 
malism. We seek self-similar solutions in which the rotational velocity and sound speed scale as 
i? -1 / 2 , where R is the radius, and consider two limiting prescriptions for the transport of angular 
momentum by convection. In one limit, the transport occurs down the angular velocity gradient, 
so convection moves angular momentum outward. In the other, the transport is down the specific 
angular momentum gradient, so convection moves angular momentum inward. We also consider 
general prescriptions which lie in between the two limits. 

When convection moves angular momentum outward, we recover the usual self-similar solution 
for ADAFs in which the mass density scales as p oc i?™ 3 / 2 . When convection moves angular 
momentum inward, the result depends on the viscosity coefficient a. If a > a crit i ~ 0.05, we 
once again find the standard ADAF solution. For a < a cr n, however, we find a non- accreting 
solution in which p oc i? -1 / 2 . We refer to this as a "convective envelope" solution or a "convection- 
dominated accretion flow." 

Two-dimensional numerical simulations of ADAFs with values of a ^ 0.03 have been reported 
by several authors. The simulated ADAFs exhibit convection. By virtue of their axisymmetry, 
convection in these simulations moves angular momentum inward, as we confirm by computing 
the Reynolds stress. The simulations give p oc i? -1 / 2 , in good agreement with the convective 
envelope solution. The i? -1 / 2 density profile is not a consequence of mass outflow. The relevance 
of these axisymmetric low-a simulations to real accretion flows is uncertain. 
Subject headings: Accretion, accretion disks — convection — hydrodynamics — turbulence 

1. Introduction 

Advection-dominated accretion flows (ADAFs, see Narayan, Mahadevan & Quataert 1998 and Kato, 
Fukue & Mineshige 1998 for reviews) were introduced to astrophysics by Ichimaru (1977) and have been 
studied extensively in the last few years (Narayan & Yi 1994, hereafter NY, Abramowicz et al. 1995, Narayan 
& Yi 1995a, 1995b, Chen et al. 1995). The recent interest in ADAFs has been stimulated to a large extent by 
the fact that these flows provide a natural explanation for many phenomena associated with low-luminosity 
accreting black holes (Narayan et al. 1998). ADAFs are also relevant for understanding black holes with 
super-Eddington accretion (Abramowicz et al. 1988). 
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NY derived an analytic self-similar solution for an ADAF which has provided considerable insight into 
the properties of these flows. A number of authors have extended this work and obtained related solutions 
(e.g. Narayan & Yi 1995a, Honma 1996, Kato & Nakamura 1998, Blandford & Begelman 1999, Manmoto 
ct al. 2000). 

NY showed that the entropy of the accreting gas in an ADAF increases towards smaller radii. They 
argued that an ADAF is therefore likely to be convectively unstable. The convective instability was confirmed 
by Igumcnshchcv, Chen & Abramowicz (1996) using numerical simulations. More detailed simulations have 
been reported by Igumcnshchcv & Abramowicz (1999) and Stone, Pringle & Begelman (1999). 

Although NY recognized the importance of convection, they incorporated its effects only schematically 
in the derivation of their self-similar solution. We present a more detailed discussion here. We show that, 
depending on whether convective turbulence moves angular momentum outward or inward, the nature of 
ADAFs may be very different. We use this insight to interpret the numerical simulations of Igumcnshchev 
and others. A related analysis is presented in the accompanying paper by Quataert & Gruzinov (2000, 
hereafter QG). 

2. Self-Similar Scalings 

Our analysis and notation closely follow the discussion given in NY. We work with height-integrated 
equations (see Narayan & Yi 1995a for an interpretation of height-integration). We assume that there 
is no significant mass outflow from the ADAF, i.e. that the mass accretion rate M is independent of 
radius R. The Kcplerian angular velocity is Qk — (GM/i? 3 ) 1 ' 2 and the corresponding linear velocity is 
vk = R£Ir = (GM/i?) 1 / 2 , where M is the mass of the accreting star (black hole) and R is the cylindrical 
radius. 

We seek a self-similar solution which satisfies the following scalings for the angular velocity ft, the 
isothermal sound speed c s and the scale height H: 

n(R) = n n K oc r- 3/2 , 

c 2 s (R) = c 2 v 2 K cx R-\ 
H(R) = c s /Q K = c R, 

where fio and Co are dimensionless constants to be determined. We write the density p as 

p(r) = p R~ a , 

so that the pressure scales as 

p(R) = pc 2 s cx R-'- a . 

The index a is equal to 3/2 in the original self-similar solution of NY; we reproduce this scaling in §4. 
However, as we show in §§5.2, 6, under appropriate conditions a very different solution is possible, which has 
a = 1/2. 

As in NY, we apply the conservation laws of mass, radial momentum, angular momentum and energy 
to solve for the various unknowns. Mass conservation requires that pvRH be independent of R, where v is 
the radial velocity. This gives the following scaling for v. 



v{R) cx R a - 2 cx R a - 3 / 2 v K . 
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In the radial momentum equation, we assume for simplicity that v 2 <C v\, which corresponds to the 
condition that the Shakura-Sunyaev viscosity coefficient a is small: a 2 <C 1. This allows us to ignore the 
ram pressure term vdv/dR. We also ignore the gradient of the turbulent pressure. The radial momentum 
equation then simplifies to a simple balance between gravity, centrifugal force and thermal pressure gradient. 
This gives 

tf K R-tfR=-±± = (a+l)<gk, 

which leads to the condition 

ng = l-(a+l)cg. (1) 
Since we require fi§ > 0, we see that Cq < 2/5 for a = 3/2 and c 2 , < 2/3 for a = 1/2. 



3. Mixing Length Convection 

We follow the model of mixing length convection developed by Grossman, Narayan & Arnett (1993, 
hereafter GNA) and use their notation, except that we replace the temperature by the isothermal sound 
speed c 2 = kTj pm p , where p, is the dimensionless molecular weight. The properties of convection depend 
sensitively on the superadiabatic gradient 

where 7 is the adiabatic index of the gas. 

To fix ideas, we consider first a non-rotating gas. The entropy gradient of the gas can be written in 
terms of AVc 2 : 

ds_ _ J_dcl _ cldp_ _ 2 
dR~ 7 -ldR pdR~ c *> AVc - 

where c p = 7/(7—1) is the specific heat at constant pressure (in units of k/pm p ). The Brunt- Vaisala 
frequency N is given by 

2_ 1 d P d ^(P 1H \ _ 5e// AY7 2 



pdRdR \ p ) c 2 

where g e ff = — (l/p)(dp/dR) is the radial effective gravity. When A^ 2 is positive, perturbations in the gas 
have an oscillatory behavior and the medium is convectively stable. However, when N 2 is negative, i.e. 
when N is imaginary, perturbations have a runaway growth, leading to convection. Convection is present 
whenever AVc 2 is positive, i.e. when the entropy decreases outward. This is the well-known Schwarzschild 
criterion. 

We assume that all mixing lengths are equal to a single length Lm, and we set the dimensionless thermal 
expansion coefficient to unity (ideal gas). We also ignore the effect of the microscopic viscosity v on the 
motion of convective eddies. This is not a very safe assumption, but it simplifies the analysis considerably. 
We believe that it will not affect the qualitative nature of the results. With these asumptions, the various 
coefficients, A, B, D, E, in GNA simplify to A = D = 0, B = E= (2/L M )- 

Borrowing results from GNA, the root mean square turbulent velocity of convective blobs is given by 

1 /2 

g e //AVcg \ _L M 2 1/2 
BEc 2 J " 2 1 ' ' 
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and the convective energy flux is 

Let us define the effective diffusion constant K c for convective energy transport by the relation 

ds 

F ' = - K ^ T dR- 

We then obtain 

K c = ^k(-N 2 f/ 2 . (2) 

We see that K c has an intuitively obvious form; it is the product of (—TV 2 ) 1 / 2 , which describes the charac- 
teristic frequency associated with convective motions (more precisely, it is the growth rate of perturbations), 
and the square of Lm, the characteristic length scale of convection. The standard treatment of mixing length 
convection in astrophysics gives a coefficient of l/4\/2 rather than 1/4 (e.g. Cox & Giuli 1968, Kippenhahn 
& Weigert 1990). Since empirical estimates of the mixing length make use of the latter coefficient, we 
incorporate this difference in our definition of the dimcnsionless mixing length Im below. 

All mixing lengths arc equal in our treatment. It is therefore natural to assume that all transport 
phenomena in the convective medium have diffusion constants of the same order as equation (2). In this 
spirit, we assume that the above expression for K c is valid also for convective angular momentum transport. 
Of course, there could be differences in the relative efficiencies of energy and angular momentum transport, 
associated for instance with different mixing lengths for the two phenomena. As we highlight in this paper, 
there are large uncertainties even with regard to the sign of the angular momentum flux. 



3.1. Convection in a Differentially- Rotating ADAF 

When there is rotation, new frequencies enter the problem and convection is no longer determined 
purely by the Brunt- Vaisala frequency TV. In the spirit of our height-integrated approach, let us consider 
displacements of gas elements in the equatorial plane. In this case, the effective frequency N e ff of convective 
blobs is given by 

N 2 eff = N 2 + K 2 , (3) 

where k is the epicyclic frequency for particle motions in the equatorial plane; for Q oc i?~ 3 / 2 , we have k = Q. 
By analogy with equation (2), we may then write the diffusion constant K c as 

K c = If {-Nt ff yi\ (4) 

QG have shown that the restriction to equatorial displacements leads to an underestimate of the 
growthrate of convective modes. In the Appendix, we present a more general analysis in which we re- 
move the restriction to equatorial displacements. In the rest of the paper, however, we work with equations 
(3) and (4). 

For a self-similar ADAF, we find that 

iV 2 = -l*A]n ( P —) = _(! + ")[!- (7 -1)<*U 2 



pdR dR 
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Ktt = {- (1 + a)[(7 + ^ (7 - 1)a] c§ + 1} n%. 

It is not immediately obvious from the above expression whether a given flow will be convectively stable or 
unstable, since that depends on the sign of which cannot be determined until the full problem is solved 
and the value of Cq is obtained. 

If N^ff < 0, we have a convectively unstable medium and the diffusion constant for energy transport is 
given by equation (4). Let us write the mixing length Lm in terms of the pressure scale height: 

L M = 2~ 1 / A l M H p , H p = -— — = — — ■ — -, 

amp (1 + a) 

where Im is the usual dimensionlcss mixing length parameter (called a in the solar and stellar convection 
literature, cf. Kim et al. 1996, Abbett et al. 1997). The additional factor of 2 -1 / 4 has been introduced 
in order to bring our formulae in line with the usual versions of mixing length convection (cf the discussion 
below eq 2). We also write the diffusion constant in the usual form 

where a c is a dimcnsionless coefficient which describes the strength of convective diffusion; this coefficient is 
similar to the usual Shakura & Sunyaev a which is used to parameterise the strength of viscosity. We then 
find that 

l 2 M f (l + tt)[( 7 + l)-(7-lM 2 A 1/2 ^ 



4V2(1 + a) 2 4 I 7 
and the convective energy flux takes the form 

c? m ds 
\Ik dR 

In the Sun, the mixing length Im is estimated to be approximately ~ 1.5 (e.g. Abbett et al. 1997). 
However, it is believed that Im should be smaller when there is a strong superadiabatic gradient (e.g. Kim 
et al. 1996); for instance, Asida (1999) estimates Im ~ 1-4 in red giant envelopes. In view of the relatively 
large superadiabatic gradient present in ADAFs, the convection in these flows is likely to be more similar to 
that found in red giant envelopes than in the Sun's convection zone. So we may expect l 2 M <~ 2 in ADAFs. 
We use this value for numerical estimates, but we note that there is some uncertainty in the estimate. 

Transport of angular momentum by convection is a complex subject and there is no consensus on how 
it operates. We consider two extreme possibilities. In one limit we assume, following NY, that convection 
behaves like normal viscosity; that is, we write the flux of angular momentum as 

• ci ^odO, 

J ° = - a <tt pR dR- (6) 

This corresponds to the assumption that the convective angular momentum flux is oriented down the angular 
velocity gradient, i.e. that convection tries to drive a system towards a state of uniform rotation, just as 
microscopic viscosity does. For Q cx i?~ 3 / 2 , it corresponds to outward transport of angular momentum. 
Equation (6) also assumes that the diffusion constant for angular momentum transport is equal to that for 
energy transport, hence the use of the same constant a c . 
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An alternative possibility is that the convective flux scales as 

• c 2 s d(nR 2 ) 

This means that the convective angular momentum flux is oriented down the specific angular momentum 
gradient, i.e. that convection tries to drive a system towards a state of uniform specific angular momentum. 
For Q oc i?~ 3//2 , it corresponds to inward transport of angular momentum. 

We concentrate on these two limiting cases in §4 and §5. However, as Kumar, Narayan & Loeb (1995) 
showed, there could be a continuum of intermediate possibilities, depending on the nature of the interactions 
between convective eddies. Therefore, we discuss in §6 the general case where the angular momentum flux 
takes the form 

Jo a c QK PU dR . (H) 

Here, the index g allows us to tunc the physics of convective angular momentum transport. When g = 1, 
we reproduce equation (6) and when g = —1/3 we reproduce equation (7). In principle, any value of g 
between these two extremes is possible. In our problem, the specific case 5 = corresponds to zero angular 
momentum transport. 



4. Outward Angular Momentum Transfer 

We assume that convection moves angular momentum outward according to equation (6) . In this case, 
the only valid self-similar solution is the one found by Narayan & Yi (1994), which has a = 3/2, p cx i?~ 3 / 2 . 
For this case, we may write the radial velocity as 

v(R) = v a v K oc iT 1/2 , 

where vo is a constant. 

We consider first the angular momentum equation. We follow Narayan & Yi (1994) and look for a 
self-similar solution in which the net flux of angular momentum vanishes. This implies that the sum of the 
angular momentum fluxes due to viscosity J v , convection J c and advection J a dv is equal to zero. Thus 

j = j v + j c + Jadv = -(a + a c )-^pR 3 ^- + pRvflR 2 = 0, 

\LK clr 

which simplifies to 

3 

vo = --(a + a c )cl. (9) 

We consider next the energy equation. Ignoring cooling, but including energy transport by convection, 
this equation takes the form 

The viscous dissipation rate is equal to (shear stress) x (rate of strain): 

0+=(a+Oe) ^wa 2 . 
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Substituting equation (1) for Q and equation (8) for w , and simplifying, this reduces to the following 
relation, 

C 2 = 9(a + a c ), (10) 

where we recall that a c is itself a function of c§ (cf. eq. 5). We note that equation (10) is equivalent to 
equation (36) in Kato & Nakamura (1998). Given the values of a and 7, this equation may be solved for Cq. 
Equations (1) and (9) then allow us to calculate the angular velocity and radial velocity, thus completing 
the solution. 

The solid lines in Fig. 1 show the variation of the convectivc coefficient a c as a function of the viscous 
coefficient a for three choices of 7: 1.6, 1.5, 1.4. We see that a c lies in the range 0.02 to 0.07, depending on 
the values of a and 7. The estimate of a c is directly proportional to l 2 M , which we have set equal to 2. 

An important result of this analysis is that we obtain a consistent solution for all values of a and for any 
7 < 5/3. In other words, the NY self-similar solution, with a = 3/2, is valid in the presence of convection 
for all values of a, provided convection transfers angular momentum outward with the same efficiency as it 
transfers energy (i.e. same a c ). This result was already obtained by NY using a somewhat simplified model 
of convection. We have presented here a more detailed analysis, using a slightly more rigorous version of 
mixing length convection. 

Note that for the NY self-similar solution the value 7 = 5/3 is special. For 7 = 5/3 and any value of 
a, we find Cq = 2/5, £!q = (no rotation) and a c = (no convection). There are no self-similar solutions of 
the NY form for 7 > 5/3. 



(5-37),„ . 45, 
^—J±(3a + a c ) + —(a + a c ) 



5. Inward Angular Momentum Transfer 

We now consider the case when angular momentum is moved inward by convection, following the 
prescription given in equation (7). To our knowledge, this case has not been considered previously in the 
theory of ADAFs. 



5.1. Self-Similar ADAF Solution 

We consider first the case when a = 3/2. The analysis proceeds very similarly as in the previous 
section. The angular momentum equation looks the same as before, except that J c is now proportional to 
—d(£lR 2 )/dR rather than —d£l/dR. Correspondingly, the vanishing of the angular momentum flux gives the 
following condition: 

v Q = --(3a- a c )cl. (11) 

In the energy equation, we note that the net viscous stress, including the effect of viscosity, is now 
proportional to (3a — a c )/2. We once again write Q + as (shear stress) x (rate of strain). Then, instead of 
equation (10), we obtain the following relation, 

(5-37)^ ^ ^ 2 _ _ ^ 5 2 



where once again a c is a function of Cq. 



(a - a c )c = (3a - a c ) 1 - -c , (12) 
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The solid lines on the right in Figure 2 show the variation of a c with a for various values of 7. As 
before, we see that a c lies in the range <~ 0.02 — 0.07. However, we now find that a consistent solution is 
available only for a greater than a certain critical a cr ni, whose value lies in the range 0.03 to 0.08. The 
value of a crit i depends on the value of 7, and is also directly proportional to l 2 M . 

When a = a crit i, the sound speed takes on its maximum value, c 2 , = 2/5 (and the rotation goes to zero, 
cf. eq 1). Knowing this fact, it is straightforward to show that 



where a c . cr i t i is the value of the convective a c when a — a cr i t i- Note once again that the solution is valid 
only for 7 < 5/3. 



What happens if a < a crit i? In this case, there is no NY- like self-similar solution with a = 3/2. A 
possible reason for the lack of solution is the following. Since advection moves angular momentum inward 
and since we have assumed (in this section) that convection also moves angular momentum inward, the only 
way to have a consistent self-similar accretion solution is for a viscosity to move an equivalent amount of 
angular momentum outward. If the parameter a is very small, the viscous flux is unable to cope with the 
inward flux due to convection, and there is no consistent accretion solution. 

However, we find that when a is small, a completely different solution, with a = 1/2, p oc i? -1 / 2 , is 
possible. This is a non-accreting solution with v = (at least in the limit of perfect self-similarity). We refer 
to it as a "convective envelope" solution, or a "convection-dominated accretion flow" (QG). 

Since v = in this solution, the advected angular momentum flux J a d v vanishes. Therefore, we must 
have J v + J c — 0, i.e. there must be a perfect balance between outward angular momentum transport via 
viscosity and inward transport via convection. This leads to the following condition on c§: 



If a is very small, then a c is also very small. In this limit, we have a flow which is marginally stable to 
convection, with c\ — 47/3(7 + 3). (This result is valid only if we restrict the analysis to equatorial motions; 
see the Appendix for a more general discussion.) 

The convective envelope solution satisfies the energy equation trivially. Since the net angular momentum 
flux vanishes, there is no shear stress, and Q + = 0. Furthermore, since v = 0, there is no advection of entropy. 
Finally, since a = 1/2, the convective energy flux F c oc R~ 2 and therefore the divergence of convective energy 
flux vanishes. This last property, which is essential in order to satisfy the energy equation, is what fixes 
a = 1/2 for the solution. 

Equation (13) allows us to solve for a c as a function of a and 7. The solid line segments on the left 
in Fig. 2 show the results. Interestingly, we find that there is a consistent solution only if a is less than a 
critical a cr i t 2- The critical value is obtained by setting Cq equal to the largest value it is allowed for a = 1/2, 
viz. 2/3 (see eq 1). This gives 




5.2. Self-Similar Convective Envelope Solution 




(13) 
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We see that a cr it2 is of the same order as a cr iti, i.e. ~ 0.05, but that the two are not exactly equal. An 
interesting feature is that 7 = 5/3 is not special for this solution. Consistent flows can be obtained even for 
larger 7. 

Note that a convection-dominated accretion flow is technically not an accretion flow at all, but a static 
configuration. In practice of course, a small amount of mass will flow into the central black hole from 
the innermost region and this will drive a small amount of accretion. However, in contrast to the ADAF 
solution described by NY, where the mass accretion rate is determined by the density and temperature on 
the outside and is not sensitive to the radius of the inner boundary, in the convective envelope the accretion 
rate is determined entirely by the conditions at the inner boundary and therefore on the location of the inner 
boundary. 

Another interesting feature is that the convective envelope solution has a significant outward flux of 
energy carried by convection. This energy flux is constant with radius and clearly originates near the center. 
Indeed, the flux originates from the (small amount of) mass which is accreted at the inner edge. A fraction 
of the binding energy of this accreted mass is diverted outward by convection and is transported to large 
radii, where it is presumably radiated in some fashion. This feature of the convective envelope solution, 
namely that energy which is generated deep in the center is transported a large distance radially before 
being radiated, is very similar to what happens in the convection zone of low-mass main-sequence stars like 
the Sun and red giants. 



6. General Angular Momentum Transfer 



For completeness, we briefly discuss the general case in which the angular momentum transfer is de- 
scribed by equation (8) with an arbitrary index g. 

If we consider an NY-like self-similar solution with a = 3/2, equations (10) and (12) are replaced by the 
more general relation 



(7 - 1) C ° + 2 C ° 9 



9g 



(3 5 - 2)(5 - 3 7 ) 2 45. 9 



'0 



(7-1) 

where a c is still given by equation (5). It is easily verified that this equation reduces to equation (10) for 
g = 1 and equation (12) for g = —1/3. We can also show that the critical ot cr i t \ becomes 

1/2 

otcnti = I - - g ) a^ cr iti = \ - - g 



20 



7 



We see that a cr u\ is positive whenever g < 2/3. This is an interesting result. It shows that the breakdown 
of the NY solution at small a is not unique to the prescription (7) for the angular momentum transfer, 
which corresponds to g = —1/3. The same thing happens even when there is no angular momentum transfer 
by convection (g = 0) or when there is angular momentum transfer outward but with a reduced efficiency 
compared to equation (6) (0 < g < 2/3). In these cases, the critical a cr n\ is smaller than when g = —1/3, 
which means that a has to be lower before the NY solution would fail. 

Similarly, when we repeat the analysis of §5.2 that led to the convective envelope solution, we find for 
a general value of g that equation (13) is replaced by 

1/2 



a = -ga c = - 



9l 2 M 
9\/2( 



3(7 + 3) ( 
47 



1 
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We see that the convective envelope solution is allowed for any g < 0, but not for g > 0. The critical a crit 2 
up to which the solution is available is 



Ot C rit2 — —gOi c ,crit2 — ^ 



9 l M ( 3 -7 X 



One interesting fact that emerges from this analysis (as well as the simpler version presented in §5) is 
that there are certain parameter regimes for which neither an NY-like self-similar solution nor a convective 
envelope solution is possible. Any flow with < g < 2/3 and a < a C riti belongs to this category. For this 
parameter range, power-law solutions, if any, would violate one of our basic assumptions. The most likely 
modification is the scaling f2 oc O^- oc i?~ 3 / 2 which we have assumed. Honma (1996), for instance, discovered 
an ADAF solution for 7 = 5/3 which has f2 oc i? -1 / 2 . This solution has been developed further by Kato & 
Nakamura (1998) and Manmoto et al. (2000). 



7. Comparison with Numerical Simulations 

In this section we compare the theory developed in the previous section with numerical results from a 
low- viscosity two-dimensional (2D) simulation of a non-radiative accretion flow. The model was calculated by 
solving the non-relativistic time-dependent Navier-Stokes equations for an accretion flow in the gravitational 
field of a point mass. All components of the viscous stress were included. The details of the numerical 
technique are described in Igumenshchev & Abramowicz (1999). 

The model considered here has a = 0.01, 7 = 5/3, inner boundary of the accretion flow at Ri n = 3R g 
and outer boundary at R out = 8x 10 3 i? g , where R g is the gravitational radius of the black hole. Mass is 
steadily injected within an equatorial torus near the outer boundary of the grid. There is no cooling in the 
accreting gas. 

Initially, there is little or no mass in the grid. As the injected mass spreads and accretes, the mass 
within the grid increases. After a period of time comparable to the viscous time scale at the outer radius, 
the accretion flow achieves a quasi-stationary behaviour, and may be considered to be in steady state. 
The "steady state" is, however, only in a time-averaged sense since the flow has convective motions which 
introduce chaotic fluctuations at any given point. We compute time- averaged properties of the flow for 
comparison with the theoretical predictions by averaging over a time equal to 44 Keplerian periods measured 
at WORg. 

Before we describe the results we would like to emphasize an important point. Stone & Balbus (1996) 
showed that if azimuthal pressure gradients are absent, then turbulence in a Keplerian disk can only transport 
angular momentum inward. By symmetry, 2D axisymmctric simulations do not have azimuthal pressure 
gradients. Because of this, each eddy preserves its angular momentum as it moves (in the absence of 
ordinary viscosity). Therefore, turbulent mixing of eddies in axisymmetric simulations tends to drive the 
system towards a state of constant specific angular momentum, and the transfer of angular momentum 
behaves like equation (7). We emphasize that this is merely a consequence of axisymmetry, and does not 
necessarily represent the properties of real convection. 

Our numerical simulations confirm the above expectation for the direction of transfer of angular mo- 
mentum. Figure 3 shows the (R, </>)-component of the Reynolds stress = (pv' R v'<) as a function of radius, 
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where the stress has been averaged over polar angle 9 as follows: 

IR^U)- j* p{Re)s[ned() ■ 

We see that, except near the boundaries (R < 5R g and R > 10 3 R g ), the Reynolds stress is negative 
everywhere in the flow. This indicates that the bulk convective motions in the gas move angular momentum 
inward. Equating the numerical estimate of ctr^ and the convective stress K c p(dflR 2 /dR)/R, wc obtain 
an estimate a c ~ ^( v R v ' r j,) / v k = 10~ 2 , where we have used H/R = 1/2. Thus, a c /a ~ 1 in the numerical 
model. 

Because the convective angular momentum transfer in the simulation follows equation (7) , we must make 
use of the analysis presented in §5 for interpreting the numerical results. Moreover, since a is very small 
(a = 0.01), we expect the simulation to reproduce the low-a convective envelope solution with p oc i? -1 / 2 
rather than the high-a NY branch of solution. To test this prediction, we show in Fig. 4 the profiles of p, 
c l/ v Ki \ v \/ v k and Q/Qk as functions of R. Each quantity has been averaged over 9 exactly as in the case 
of the Reynolds stress. 

The simulation gives a density profile with p oc i?~ 1//2 , as expected from the theory. A similar R^ 1 / 2 
variation was seen also in low-viscosity ADAF simulations by Igumenshchev & Abramowicz (1999) and 
by Stone et al. (1999). The latter authors considered a variety of viscosity prescriptions, of which Run 
K corresponded to the usual a prescription; this run exhibited p oc R~ x / 2 behavior, as expected for the 
convective envelope. It may be worth emphasizing that the R" 1 ^ 2 density profile in the convective envelope 
solution (which is very different from the i?~ 3 / 2 profile of the self-similar NY solution) is not a consequence 
of mass loss. In the analytic work presented in this paper, for instance, we explicitly assume that there is 
no mass loss. Both Igumenshchev & Abramowicz (1999) and Stone et al. (1999) found that for low values 
of a there were no powerful unbound outflows in their numerical simulations. 

In this context, note that because the medium is convectively turbulent and is not accreting, roughly 
half the mass at any radius at any given time will be flowing in and half will be flowing out. We could 
estimate a "mass inflow rate" Mi n (R) by adding up all the inflowing gas elements. This would give Mm <~ 
4npaRH/2 where a is the root mean square velocity of the turbulent eddies. In §3 we showed that a scales 
as L M (-^V 2 // ) 1/2 Rr 1 ! 2 . Therefore, for p oc i?" 1 / 2 and H oc R, wc find that M in oc R. It might be 

tempting to identify the rapid decline of M in with decreasing radius as evidence for a massive outflow. But 
there is no reason to think that the outflowing gas at any radius in the interior of the flow will escape and 
flow out to infinity. It is more than likely that the outflowing gas forms part of convective eddies, so that 
the mass that is flowing out at a certain time will later turn around and flow in. If such is the case, then the 
quantity Mi n is not very useful. 

The analytic convective envelope solution predicts that M = and v — 0. As Fig. 4 shows, the mean 
radial velocity in the simulation is extremely small, but it is not exactly zero. The small non-zero value is 
because the simulated flow has an inner boundary where some mass can flow into the black hole. Thus there 
is a small M which leads to a finite v. The scalings discussed in §2 show that for p oc R^ 1 / 2 , we expect the 
velocity to scale as v oc vk/R- This is confirmed in Fig. 4. 

The profiles of c 2 {R) and Q(R) in Fig. 4 show the predicted power-law behaviors of these quantities for 
10i? 9 < R < 10 3 R g . However, a comparison of actual numerical values does not give good agreement. From 
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equations (1) and (13), for a = 1/2, 7 = 5/3 and very small a, the theory predicts 

c 2 47 ft /3- 7 \ 1/2 

^ = ^7^=0.48, _=(-_I) =0.53. 



^ 3(3 + 7) ' ' Ojf \3| 7 

The numerical simulations give c 2 s /v 2 K = 0.26 and = 0.9. As QG have shown, this discrepancy is the 

result of the height-integrated approximation. We discuss in the Appendix an improved theory that avoids 
this approximation. The theory predicts (? s lv 2 K = 0.29 and Q/Qk = 0.75, in better agreement with the 
simulations. 

An important feature of the convective envelope solution is that there is a strong outward flux of energy 
F c transported by convection. The corresponding luminosity is equal to E c = 4ttRHF c , which is independent 
of radius. Comparison of the theoretical estimate with numerical results from the simulation show quite good 
agreement. Substituting the numerical value of pR}l 2 at R = 100i? g in the analytical expression for F c , 
we estimate E c ~ 0.2a c in units of Mc 2 , where M is the mass accretion rate into the black hole. On the 
other hand, the numerical simulations give E c — E tot — E v ~ 2 x 10~ 3 , where E tot is the total outward 
energy flux and E v is the energy flux due to viscosity. Equating the two estimates of the convective flux we 
get a c ~ 10~ 2 . This estimate of a c agrees with our previous estimate from the Reynolds stress, which is 
encouraging considering the uncertainties (physical and numerical) of the two estimates. Thus, we feel that 
the convective envelope model is a good description of numerically simulated 2D ADAFs at low a. 



8. Discussion 

The main result of this paper is that, in addition to the NY self-similar ADAF solution which has density 
p oc i?~ 3 / 2 , there is a second power- law solution, a radiatively inefficient "convective envelope" solution or 
"convection-dominated accretion flow," in which p oc R~ x l 2 . The second solution corresponds to a static 
envelope in which the mass accretion rate is very small. Indeed, M vanishes in the limit when the inner edge 
of the accretion flow moves down to 0. The solution has well-developed convection and is perfectly balanced 
in the sense that it transfers as much angular momentum outward via viscosity as it transports inward via 
convection. There is, however, a net outward flux of energy. 

The convective envelope solution is possible only when two conditions are simultaneously met: 

(i) The viscosity coefficient a must be fairly small. For a reasonable choice of the mixing length parameter 
(I'm ~ 2), we find that we require a < a cr i t 2 ~ 0.05 (see Fig. 2). The improved analysis presented in the 
Appendix modifies the value upward. 

(ii) Convection must transport angular momentum inward (g = —1/3, §5) or not transport angular momen- 
tum at all (g — 0, §6). If it transports angular momentum outward, it must do so with less efficiency than 
it transports energy; specifically, we require g < 2/3 (§6). 

The numerical simulations discussed in §7 and in Igumenshchev & Abramowicz (1999) and Stone ct 
al. (1999) satisfy both of the above conditions. In particular, condition (ii) is automatically satisfied by all 
numerical simulations of ADAFs carried out so far since all these simulations have been axisymmetric. The 
axisymmetry guarantees that there are no azimuthal pressure gradients and so convective eddies transport 
angular momentum inward (Stone & Balbus 1996). We find that the numerical results agree quite well with 
the analytical results (Figs. 3, 4, §7). 

In all cases where either condition 1 or 2 is violated, the only self-similar solution we have been able to 
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find is the solution derived by NY, in which p cx R 3 / 2 . For given outer boundary conditions, this solution 
leads to a considerably larger mass accretion rate than that obtained with the convective envelope solution. 

Igumcnshchev & Abramowicz (1999) found that convection plays a negligible role when a is <; 0.1. 
This is consistent with our analysis, which shows that the convective a c is much less than a in this regime 
(Figs. 1,2). Indeed, our analysis overestimates a c for such models because the mixing length formalism 
we employ is a local steady state theory which assumes that convection has achieved saturation amplitude. 
However, the convective turnover time, which is given by t conv — (— ./V 2 ^) -1 / 2 is of order l/a c flx, whereas 
the accretion time in the NY solution is of order t acc <~ 1/aQ.x ■ Thus t conv /t acc ~ a/a c , and this becomes 
large as a increases. The net effect would be to reduce a c below the values shown in Figs. 1 and 2. Thus, we 
expect convective motions to be negligibly small for a <; 0.1, as confirmed by Igumcnshchev & Abramowicz 
(1999). (Stone et al. 1999 did not explore models with such large values of a.) 

The range a ~ 0.1 — 0.3 is found to be very interesting in Igumenshchev & Abramowicz's (1999) 
simulations. For a = 0.3, 7 = 1.5, they found that there is a pure inflow, while for other choices, e.g. 
a = 1, 7 = 1.5, they obtained a stable bipolar outflow. Igumcnshchev (2000) studied a model with a = 0.1, 
7 = 5/3 and found a global meridional circulation pattern accompanied by a surprising unipolar supersonic 
outflow. The flow has some resemblance to the convective envelope solution discussed in this paper. Only 
a small fraction of the circulating matter is directly accreted by the black hole, and the energy required to 
support the circulation is extracted from the infalling mass with an efficiency <~ 10~ 3 — 10~ 2 . The circulation 
transports the energy to large radii, where the author argues it could be radiated via bremsstrahlung. It 
would appear that in this model, viscosity is large enough to suppress convection on scales ^ R, but global 
circulation on a scale ~ R still survives, and this transports energy. 

The parallel between this solution and our convective envelope solution is fairly strong. The convective 
envelope again converts a fraction of the binding energy of the accreting mass into an outward flux of energy, 
which is carried by convection. This energy somehow has to be got rid of at the outer boundary. One way 
in which the energy could be eliminated from the system is by driving a slow (but perhaps massive) outflow 
from the outer boundary. Another is that the energy could be radiated by gas on the outside, as suggested 
by Igumenshchev (2000). 

The latter is quite plausible. Since p cx i? -1 / 2 and T oc J? -1 , the bremsstrahlung cooling rate per unit 
volume varies as Q^ rems cx p 2 T x / 2 cx i?~ 3 / 2 . The cooling time scale thus goes as t coo i cx pT/Q^ rems cx R°. In 
comparison, since v cx i?~ 3 / 2 , the accretion time scale goes as t acc — R/v cx R 5 / 2 , and the ratio t acc /t coo i cx 
i? 5 / 2 increases very rapidly with radius. Therefore, one naturally expects a highly advection-dominated 
convection zone (presumably well described by our convective envelope solution) sandwiched between an 
inner energy producing zone close to the black hole and an outer radiating zone. The radiating zone would 
have very specific spectral signatures which would be interesting to explore. 

Is the convective envelope solution relevant for real accretion flows? This depends on the answers to 
two questions. 

The first question is: What is the value of a in real ADAFs? If a > a crit 2 ~ 0.05, then the convective 
envelope solution is not possible and the second question below is irrelevant. Empirical estimates of a from 
fitting the spectra of black hole systems suggest that a is large. Narayan (1996) obtained a <~ 1, which Esin, 
McClintock & Narayan (1997) later revised toa~ 0.2 — 0.3. Quataert & Narayan (1999) considered models 
with a = 0.1, but they did not verify that such low values of a are consistent with the constraints that 
Narayan (1996) and Esin et al. (1997) considered, specifically the existence of luminous black hole X-ray 
binaries in the "low" spectral state. More work is needed before one can state with confidence whether or 
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not values of a & a cr i t 2 are consistent with observations. 

On the theoretical front, although there is good reason to believe that "viscosity" in differentially- 
rotating accretion flows is produced by magnetic stresses generated by the Balbus & Hawley (1992, 1998) 
instability, the theory has not developed to the point where the value of a can be estimated. Numerical 
MHD simulations give values over a wide range, a ^ 0.01 toa~ 0.6 (Balbus & Hawley 1998). Furthermore, 
all simulations so far have been done on thin accretion disks, and it is unclear whether those results are valid 
for the much thicker ADAFs. Numerical MHD simulations of accretion flows under ADAF-like conditions 
would be very worthwhile. The simulations will need to be done in three dimensions rather than two, since 
even in the thin disk case, there are large differences between 2D and 3D simulations (Balbus & Hawley 
1998). 

The second question is: Does convection in differentially-rotating ADAFs move angular momentum 
outward or inward? 

Ryu & Goodman (1992) showed that linear modes in a convectively unstable thin accretion disk transfer 
angular momentum inward. Stone & Balbus (1996) used numerical simulations to study the non-linear 
version of the problem. They found that angular momentum was either transported very weakly inward 
or not at all. These are significant results, but their relevance to ADAFs is a little uncertain. As NY 
emphasized, the entropy gradient in a thin disk model is in the vertical direction whereas the angular 
momentum gradient is in the horizontal direction. In ADAFs on the other hand, both the entropy gradient 
and the angular velocity/ angular momentum gradient are in the radial direction. This might conceivably 
cause some differences in the physics. 

Kumar, Narayan & Loeb (1995) analysed convective angular momentum transport in a differentially 
rotating medium using a mixing length formalism. They found that, while small amplitude linear perturba- 
tions transport angular momentum down the specific angular momentum gradient, in agreement with the 
result of Ryu & Goodman (1992), nonlinear saturated convection generally behaves very differently. Indeed, 
in the nonlinear regime, the angular momentum transport depends critically on the nature of the interactions 
between convective eddies, which Kumar ct al. modeled via a scattering term. By tuning their scattering 
function, they could reproduce a wide variety of behaviors, covering the entire range from g = — 1/3 to 5 > 1 
(in the language of the present paper). 

Stone & Balbus (1996) analyzed the basic hydrodynamic equations in a turbulent Keplerian disk and 
showed that, when azimuthal pressure perturbations are small, there can be no net outward transport of 
angular momentum. As already mentioned, this powerful result is the reason why axisymmetric numerical 
simulations of convecting disks always move angular momentum inward. It also throws considerable light 
on the Kumar et al. (1995) work since the eddy scattering function invoked by these authors implicitly 
involves azimuthal pressure fluctuations. It is the scattering that enables their model to transport angular 
momentum outward. 

Despite these theoretical studies, the basic question of how angular momentum transport actually oper- 
ates in a convection-dominated accretion flow remains very much open. Specifically, it is unclear how impor- 
tant azimuthal pressure perturbations are in these systems, and whether angular momentum is transported 
outward or inward. We see only one way to answer this question: full three-dimensional numerical simula- 
tions. If 3D simulations are too expensive, one should at least carry out nonaxisymmetric two-dimensional 
simulations in cylindrical geometry. 

One empirical fact is worth noting; the convection zone in the Sun is closer to being in a state of constant 



- 15 - 



angular velocity than a state of constant angular momentum. Indeed, in the equatorial plane, the angular 
velocity actually increases with increasing radius, which motivated Kumar et al. (1995) to consider models 
with g > 1. One of the two conditions, (i) and (ii), mentioned at the beginning of this section must be 
violated in the Sun. One possibility is that the effective viscous a in the Sun is large enough to counter 
any inward transport of angular momentum by convection. This is reasonable since the convective motions 
in the Sun are not very energetic, and a does not have to be very large to counter the effect of convection. 
Alternatively, perhaps convection in the Sun behaves as in equation (6) and drives the system to near-zero 
dtt I dR. The Sun is certainly not a good comparison for an ADAF — a more rapidly rotating convective star 
would be closer. Nevertheless, the fact that the Sun has achieved an equilibrium configuration that is very 
different from a constant specific angular momentum state must be treated as a possible clue as we grope 
towards a better understanding of advection-dominated flows. 

In closing, we note that the key point of this paper is that there is considerable uncertainty associated 
with angular momentum transport by convection, and that this leads to two very distinct configurations for 
an ADAF. There is a related, though less uncertain, issue connected to energy transport. We considered in 
this paper convective transport, where energy is transported down the entropy gradient. In addition, there 
could also be conduction which transports energy down the temperature gradient. Conduction has been 
studied by Gruzinov (1999) in the context of Bondi accretion. It would be of interest to extend Gruzinov's 
analysis to ADAFs and to investigate models in which convection and conduction are both present in an 
accretion flow. 
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Appendix A 

The model of convection which we employed in the main text of the paper is somewhat limited because of our use of 
height-integrated equations. In particular, when we calculated the diffusion constant associated with convection, we set the 
characteristic growthrate of convective motions equal to {—N 2 ^) 1 / 2 , where N^jf = N 2 + re 2 , and N is the Brunt-Vaisala 
frequency and re = U is the epicyclic frequency (see eqs 3 and 4). This is appropriate for a one-dimensional flow in which 
convective blobs experience displacements only in the equatorial plane. However, as QG have pointed out and as is well-known 
for rotating fluids (e.g. Tassoul 1978, Bcgelman & Meier 1982), in the full 2D flow the most unstable modes are not in the 
equatorial plane. In this Appendix we present a modified analysis in which we identify the most unstable mode and use its 
properties to estimate the effects of convection. We draw heavily on the work of QG (specifically, the Appendix of their paper). 

We consider an axisymmetric self-similar 2D flow in polar coordinates: rd. We assume that the density varies as r~ a (QG use 
n instead of a) and that the azimuthal velocity and the sound speed c s vary as r -1 / 2 . Using QG's results as a guide, we 
assume that v$ and c s at a given r are both independent of 8. (QG showed this only for the convectively marginally stable 
system, but we assume that the result is valid also for the more general situation we consider.) In terms of the Keplcrian 
velocity vjf = (GM/r) 1 / 2 , we may write 

v 4 ,(r,6) = Sl Q v K {r). (Al) 

Consider a blob of material at some arbitrary rd and imagine imposing a unit displacement: dr = cosx, rd0 = sinx- 
Following the analysis presented by QG, the effective frequency N e ff associated with the dynamics of the blob is given by 



n 2 



A cos x + B sin X + C cos x sm X: (A2) 



where the quantities A, B and C are (we have translated the expressions given in QG to the notation of this paper): 

A = 1 - a( 7 - 1) + [a( 7 - 1) - (7 + l)]n§, (A3) 

(a + 1)(7 - i)cot 2 eni , , , s 

B= (1-n 2 ) 27 °' ( } 

C = 2[a(7- 1) - (7 + 1)] cot 0Q 2 ,. (A5) 
QG show that there are additional terms in B and C which depend on dv<p /d6, but these terms vanish under the assumption 
that v$ is independent of 8. 

In the height-integrated problem considered in the main text, there is a unique frequency associated with convective motions, 
namely N eff = (N 2 + re 2 ) 1 / 2 . In the more general problem considered here, there is a multitude of frequencies; indeed, the 
frequency is a function of two angles, the polar angle 8 corresponding to the location of a blob, and the angle x associated with 
the direction of displacement of the blob. (The analysis in the main paper corresponds to the specific choice 8 = 7r/2, \ = 0.) 

In the two-dimensional space of 9 and x, the value of — -^e// var i cs as a function of the two angles. If —N 2 ^ < (i.e. 
7V 2 j,y > 0) for all choices of the angles, then the medium is convectively stable. In this case, the convective diffusion constant 
K c is zero and so is the dimenionless constant a c . We are not interested in such flows. Rather, we are interested in those flows 
in which —N 2 ^ > for at least some choices of 8 and x- If we determine the largest positive value of — N 2 ^ in such a flow, 

then (-^e//) 172 would 

give the growthrate of the most unstable convective mode. We assume that this most unstable mode 
dominates the convection. 

To determine the largest positive value of —N 2 ^, we first differentiate (A2) with respect to x and set the result to zero. We 
find that the most unstable mode at a given 8 corresponds to displacements that satisfy 

C ° S2X= [{A ~B) 2 + C 2 ]i/ 2 > Sia2 ^= [{A . B) 2 + c2]i/ 2 - (^ 6 ) 
Substituting this value of \ back in (A2) we obtain the maximum value of (—N 2 ^) at each 8: 

(-^T 1 ) =\((A + B) + [(A-B) 2 +C 2 ] 1 ' 2 ). {AT) 

V K J o,max 

Maximizing this quantity next with respect to 8, we obtain the overall maximum value of (~N 2 ^). We find that for all flows 
that are convectively unstable, i.e. which have —N 2 ^ > over at least some range of 8 and x, the maximum positive value of 
—N 2 jj occurs at 8 = 0. Now, as 8 — * 0, A remains finite, but B and C both diverge, with B diverging faster than C (cot 2 8 
versus cot 8). This simplifies matters, and so we find that 

• 

K / max \ / 9^0 
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Both in (A7) and (A8), the corresponding mode is convectively unstable only if the quantity on the right hand side is positive. 
We focus only on such cases. 

As noted by QG, the simple "equatorial" analysis presented in the main text of this paper corresponds to using just A in 
equation (A8). The more complete analysis presented here shows that the correct expression corresponds to the quantity 
(A — C 2 /4B). This is larger than A because B is negative. Thus, for instance, it is possible for A to be negative — which 
means that the medium is stable to equatorial displacements — and yet displacements at non-equatorial values of 9 could be 
unstable, so that the medium as a whole could be convectively unstable. To evaluate (A8), we note that although B and C 
diverge as 9 — » 0, the ratio C 2 /B does not, and so (—N 2 ^) ma x remains finite as 9 — > 0. The full expression is 

\ =i_o( 7 -l) + [o( 7 -l)-(7+l)]ng 



[q( 7 -l)-( 7 + l)] 2 ng(l-^) 

27(l-ng) + (o + l)( 7 -l)ng ■ 1 ' 

Having determined (—N 2 ^) max , the rest of the analysis proceeds exactly as in the main paper, except that we replace 
equation (4) for the convective diffusion constant by 

K c = ^k(-N 2 ff ) 1/2 . (AW) 

4 V C J J / max 

Note that this result is valid only if {—N 2 ^) max > 0; otherwise, K c = 0. Just as equation (4) is replaced by (A10), equation 
(5) is also correspondingly modified. None of the other equations needs to be modified. The radial momentum equation, for 
instance, is still given by equation (1). For NY-likc self-similar flows with a = 3/2 we make use of equation (10), while for the 
convective envelope solution with a = 1/2 we use equation (12). 

The dashed lines in Figs. 1 and 2 show some numerical results. We see that the results are qualitatively similar to what we 
obtained with the simpler theory presented in the main text (shown by solid lines in the plots). The value of ct c is larger now 
compared to the height-integrated theory, as one might expect. The critical a cr it2 up to which the convective envelope 
solution is possible is also larger by a few tens of per cent. Interestingly, the other critical a cr it\ is unaffected. When 
a = a cr ui, the solution has no rotation and the flow is perfectly spherical. In this limit, equatorial displacements capture the 
full story and the more detailed analysis presented here adds nothing. 

The problem considered by QG corresponds to a convectively marginally stable accretion flow in which a — » 0, a c ~ * 0. They 
show that in this limit the Bernoulli parameter, 

Be=\vl + \v 2 + -L-c 2 s -v 2 K , (All) 

is zero. We confirm this result. The reader is referred to NY, Narayan &: Yi (1995a) and Blandford & Begelman (1999) for a 
discussion of the significance of the Bernoulli parameter. For non- vanishing a and a c , we find that the Bernoulli parameter is 
always positive. For a = a cr u2, for instance, we obtain Be ~ 0.3i>^-. 

The positivity of the Bernoulli parameter for finite a c in the convective envelope solution is fairly easy to understand once we 
realize certain properties of these flows. First, in these flows v = 0, and and c s are independent of 9 (according to the 
analysis presented by QG which we have borrowed). Therefore, even though the flow rotates and is non-spherical, Be is 
nevertheless independent of 9 and is a function only of r. Second, the flows we are interested in are convectively unstable. In 
order to have a convective instability, the Bernoulli parameter has to be a decreasing function of increasing radius. Finally, 
the flows are self-similar, which implies that Be must be proportional to The latter two conditions can both be 

satisfied only if Be is positive. Therefore, for finite ct c , we expect a positive value of Be. In the limit a c — » considered by 
QG, the flow is marginally stable to convection. This means that Be must be independent of r, which can be reconciled with 
self-similarity only if Be = (as QG found). 
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Fig. 1 . — Variation of the convective coefficient a c as a function of the viscosity coefficient a for three values 
of the adiabatic index 7. The calculations correspond to the case when convection moves angular momentum 
outward according to equation (6). The solutions are of the self-similar form derived by NY, in which the 
density scales as p oc i?~ 3 / 2 . Note that there is a consistent solution for all values of a. The solid lines 
correspond to the height-integrated version of the theory described in the main paper, and the dashed lines 
correspond to the analysis described in the Appendix. 
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Fig. 2. — Variation of the convective coefficient a c as a function of the viscosity coefficient a for three 
values of the adiabatic index 7. The calculations correspond to the case when convection moves angular 
momentum inward according to equation (7). The lines on the right refer to a solution of the self-similar 
form derived by NY, in which the density scales as p oc i?~ 3 / 2 . This solution is only available for a greater 
than a critical value a cr iti, which depends on 7. The lines on the left correspond to the convective envelope 
solution discussed in §5.2 in which p oc i? -1 / 2 . This solution is only available for a less than a critical value 
otcrit2i which again depends on 7. The solid lines correspond to the height-integrated version of the theory 
described in the main paper, and the dashed lines correspond to the analysis described in the Appendix. 
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Fig. 3. — Shows the (R, </>)-component of the Reynolds stress tensor, <TR<p 7 as a function of radius in an 
axisymmetric numerical simulation of an ADAF with a — 0.01 and 7 = 5/3. The Reynolds stress has been 
averaged over the polar angle 9 (as explained in the text), and normalized to v\. Apart from two regions 
near the boundaries, R < 5R g and R > 5 x 10 3 i? g , the stress is negative. This confirms that convection in 
axisymmetric simulations moves angular momentum inward, as argued by Stone & Balbus (1996). 
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Fig. 4. — Profiles of p (solid line, upper left panel), c 2 s /v 2 K (upper right), |w|/w_ff (solid line, lower left) 
and VL/VLk (lower right), shown as functions of radius, for an axisymmetric convectively unstable accretion 
simulation with a = 0.01 and 7 = 5/3. All quantities have been averaged over polar angle 9 and time. 
Except near the boundaries (R < 10R g and R > 10 3 i? g ), the profiles of p, (? s jv 2 K and show the 

power-law behaviors predicted by the self-similar convective envelope solution; also the radial velocity v is 
small as predicted. In the upper left panel the dashed line corresponds to the analytical scaling p oc i? -1 / 2 . 
Similarly, in the lower left panel the dashed line corresponds to the predicted scaling \v\/vk oc R~ x . 



